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Abstract 

A theory of vibrational energy relaxation based on a semiclassical treatment 
of the quantum master equation is presented. Using new results on the semi- 
classical expansion of dipole matrix elements, we show that in the classical 
limit the master equation reduces to the Zwanzig energy diffusion equation. 
The leading quantum corrections are determined and discussed for the har- 
monic and Morse potentials. 



I. INTRODUCTION 



Vibrational energy relaxation (VER) is of fundamental importance for many chemical 
reactions in solution^!. VER studies in molecules like I 2 , Hgl, 2 , and N 2 in various solvents 
by means of pulsed laser techniques cover a range of relaxation times from picoseconds to 
secondJi. A detailed theoretical description is often complicated by the fact that various 
energy relaxation processes contribute. While in the high frequency regime a quantum 
mechanical treatment, e. g., on the basis of the isolated binary collision model is natural, the 
low frequency regime has for the most part been treated in terms of classical theories such as 
the Landau-Teller model, the generalized Langevin equation, or by molecular dynamicsH. 
These approaches are justified for high temperatures and heavier molecules. 

However, in some of the above mentioned molecules in solution as well as for adsorbates 
on surfaces, molecules in cryogenic matrices, or molecular crystals, purely classical theories 
are not sufficient and quantum effects cannot be disregarded. This is also due to the fact 
that experiments have to be performed at low temperatures where only few vibrational levels 
are occupied. For these systems, VER was mostly treated on the basis of Fermi's golden 
rule a„d the quantum master equatio M. 

One of the questions arising in this context is the connection between the quantum 
mechanical theory of VER and classical energy relaxation. To address this problem, we 
give here a semiclassical analysis of the quantum master equation. The Zwanzig energy 
diffusion equation! is recovered in the classical limit and the dominant quantum corrections 
are determined explicitly. 

The paper is organized as follows. In section II we introduce the microscopic model which 
is a nonlinear oscillator coupled to a bath of harmonic oscillators. This model is known to 
lead in the classical limit to the generalized Langevin equation. Molecular dynamics simu- 
lation suggest that the model can also be used for rather complicated relaxation processes, 
such as e.g., a binary molecule in a Lennard- Jones fluid, if the parameters of the harmonic 
oscillator bath are chosen appropriately^. We briefly discuss the classical limit of the model, 
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and its weak-coupling limit, which leads to the Zwanzig energy diffusion equation. Further, 
the fully quantum mechanical treatment of VER on the basis of a master equation is given. 

Section III is devoted to the semiclassical expansion of the master equation. This re- 
quires, apart from known expansions of wave functions and energy levels, explicit results on 
semiclassical dipole matrix elements. On the basis of these results we then show that the 
master equation reduces to the classical energy diffusion equation in the limit H — ► 0, and 
determine the leading quantum corrections. Section IV contains explicit results for the cases 
of the harmonic and Morse potentials, and in section V we present our conclusions. 



II. THE MODEL AND FORMULATION OF THE PROBLEM 

Vibrational motion of molecules can be modeled as the motion of a particle in an anhar- 
monic potential well V(q). To account for dissipation we couple the particle to a bath of 
harmonic oscillatorsll. Of course, the harmonic oscillator bath is a crude simplification of the 
real interaction of molecular vibrations with the solvent. However, inasmuch as the back- 
action of the solvent can be treated within linear response theory, the bath can be modeled 
in terms of effective bath oscillators. In fact, this simple model describes the classical dy- 
namics of a binary molecule in a Lennard- Jones fluid quite accurately!. For the problem of 
energy relaxation we expect the harmonic oscillator bath to be a reasonable starting point. 
With these premises, the Hamiltonian H is composed of the oscillator part 

B. = ^ + V( q) , (1) 



the bath Hamiltonian 



N 2 

tt \ ~* Pn , m n 2 2 /rA 

^b = L^- + ~iru n x n , (2) 

n=l z '" l n z 



and an interaction part 



N N 2 



tfsB = -gE^ + g 2 £^- (3) 

n=l 71=1 Z " ( « U 'n 
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It is well knownll that the associated classical equation of motion is the generalized Langevin 
equation 

Mq{t) + V'(q(t)) + M f ds 7 (t - s)q(s) = £(t) (4) 

Jo 

where the prime denotes differentiation with respect to the argument, and 

is the damping kernel. Finally, £(t) is a stochastic force having the properties 

<*(*)> = (6) 

and 

mm) = MksT^t). (?) 

We will be interested in the motion at weak damping. When uq is the frequency of small 
undamped oscillations about the minimum of the potential well, a typical damping strength 
is given by 



oc 



7 C = / dt^(t)cos(u t). (8) 





Further, we may introduce a typical memory delay time 

J °°dtt-f(t) cos(uj t) 
Tc J3»(ft7(t)cos(wb0' 

The region of damping parameters examined in this article is then characterized by 

7c < u , l/r c . (10) 

Let us briefly recall the classical treatment of this weak damping limit. Under the con- 
ditions stated, the unperturbed energy 

E=Mq 2 + V(q) (11) 



of the particle will be an almost conserved quantity, and it is advantageous to rewrite the 
classical Langevin equation in terms of energy and phase variables. 
The unperturbed trajectory q(E, t) of energy E satisfying 

Mq + V\q) = (12) 

may be written as a Fourier series 

oo 

q(E,t) = £ Q,(£)e^ (t) (13) 

Z = — CO 

where 

= uj(E)t. (14) 

Here 2tt/uj(E) is the oscillation period of a trajectory of energy E in the potential well. The 
phase will be chosen such that for $ = the trajectory is at a turning point. Then the 
Fourier coefficients are real and satisfy 

Q.t{E) = Qi{E). (15) 

Now, relation flT3| ) and 

CO 

q(E,t)= J2 Hu(E)Qi(E)e m ® (16) 

l = — CO 

define a mapping from the variables q, q to the variables E, $. By means of this transfor- 
mation, the Langevin equation may be written in terms of energy and phaseB§0. These 
transformed equations may be solved perturbatively to second order in the coupling to the 
bath. Within the time window 

l/u ,r c <t < l/ 7c (17) 

the solution can be matched with the solution of a Fokker-Planck equation for the energy 
and phase distribution function^. The structure of the resulting Fokker-Planck process 
allows one to average over the phase variable yielding the one- dimensional Fokker-Planck 
equation 



p W ^a (B )(i + * b tA)^p (b) . (18) 

where 

oo 

A(E) = 2ttMuj(E) l 2 Ql{E) x i(iluj(E)) (19) 

l=— oo 

is the energy relaxation coefficient. Here 

POO 

j( z ) = j dte- zt -f(t) (20) 

is the Laplace transform of the damping kernel. Eq. fll8|) is the energy diffusion equation 
first derived by ZwanzigS. 

We now turn to the problem of energy relaxation in a weakly damped quantum system. 
Introducing the energy eigenstates of the unperturbed system 

H s \n) = E n \n), n = 0,1,2,.... (21) 

we start from the Pauli master equation 

oo 

Pn{t) = Wn,mPm{t) ~ W m>n p n {t)] (22) 
m=0 

where p n {t) is the probability to find the system in state n at time t. The transition rate 
from state m to n may be calculated by Fermi's golden rule 

1 r°° 

W n , m = -2 / dtTr B [(n\H SB (t)\m)(m\H SB \n)p e £] . (23) 
n J-oo 

Here Tr B denotes the trace over the bath states, p B q = e _/3HB /TrB[e~^ B ] is the equilibrium 
bath density matrix, and H SB (t) = e ^ Hii+HB)t/h H SB e~^ Hs+HB)t/h . Using Eq. ©, we obtain 

N 

.,-(n\a\m) 2 



1 N fOO 

W n . m = -2 (n\q\m) 2 ^ c k c x \ dte i ^ En ~ Em ^ t ^ h Tr B [x k (t)x lP c B q ] . (24) 
n k,i=i J -°° 



In the weak damping limit considered here, the so-called counter-term, that is the second 
term of H$ B in Eq. @, may be disregarded. The master equation (|22| ) is appropriate for 
weak damping satisfying condition (|l~0|) . In addition, for the quantum system a further 
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frequency scale set by the lowest Matsubara frequency v = 2nkBT/h becomes relevant. As 
can be seen from the derivation of the Pauli master equation, bath correlation functions are 



assumed to decay on a time scale much shorter than l/j c (see for instance Ref. |TT|). For 
low temperatures where v < l/r c bath correlations decay on the time scale 1/v. Hence, in 
addition to Eq. (jnj) we have to require 

7c < v (25) 

which is, however, usually obeyed for systems of relevance in chemical physics. 

Due to the linearity of the bath, the transition rate may be evaluated exactly yielding 

W n<m = ^(n\q\m) 2 D(E n - E m ). (26) 

Here 

D(E) = 2MEn(E)Re A /(i\E\/h) (27) 

is an effective, occupation weighted spectral density of the bath with the Bose occupation 
function 

"(*0 = i^rrr ( 28 ) 

Note that the coupling to the bath is completely described in terms of the Laplace transform 
of the damping kernel. Since the detailed balance condition 

W m ,„ = e^ En - E ^W n , m (29) 

is clearly obeyed, it is readily seen that Eq. (p2]) has the stationary solution 

p n = \e' &En (30) 



where 



Z = e~ PEn (31) 

71=0 



is the quantum mechanical partition function. 



III. SEMICLASSICAL EXPANSION OF THE PAULI MASTER EQUATION 



We now want to regain the Zwanzig energy diffusion equation from the Pauli master 
equation and determine the leading quantum corrections. The results of this somewhat 
lengthy calculation are summarized at the beginning of section IV. We start by introducing 
the probability distribution function 



P(E) := £ 6(E - E n ) Pn (32) 

n=0 

by means of which expectation values are simply calculated as 

/oo 
dEP(E)f(E) (33) 
-oo 

for arbitrary f(E). In terms of the function N(E) which counts the number of energy levels 
below energy E 

oo 

N(E) = J2®(E-E n ), (34) 

n=0 

the probability P(E) can be written 

P(E) = p(E)p N{E) (35) 

where 

P(E) = ^± = ±6(E - E n ) (36) 



is the level density. Now, combining Eqs. fl22|) and (j32|) we find 

oo 

P(E)= £ [S{E - E m ) - 6{E - E n )]W m>nPn 

71,771=0 

OO oo 

= E E i S ( E - E n+i) - S(E - E n ))W n+l>nPn . (37) 

n=0 l=—n 

This equation will be evaluated in the semiclassical limit in the sequel. 
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A. Semiclassical energies 



To proceed we need to determine the semiclassical energy eigenvalues and dipole matrix 
elements. We start from the semiclassical (WKB) expansion of the wave function 



where C(E) is a normalization factor. Inserting this ansatz into the Schrodinger equation 
and expanding in powers of h, one arrives at the usual eikonal expansion 

S(E, q) = S (E, q) - ihS^E, q) - h 2 S 2 (E, q) + ih 3 S 3 (E, q) + . . . (39) 

where the Si(E, q) can be calculated recursively with the resultll! 

^-S (E,q)=p(E,q) = pM[E - V(q)], 
| 5l(£ ,,) = -I^l„ p(B , 5 ), 

~l 2 



— S 3 (E q) = 1 - 1 - S 2 (E q) 
dq ' 2dqp(E,q)dq 

To second order in H the wave function is then given by 



§- q P(E,q) 



(40) 



H 2 dS 2 (E,q) 



0{% 6 )\. (41) 



2p(E, q) dq 

The solution to the two turning point problem we are facing here is a linear combination 
of these right and left traveling waves. Details of the analytic structure of these functions 
are discussed in the literature@. Here we only mention that the quantization condition is 
obtained as 

1 ldq dS ^' q) =2mn, n = 0, 1,2, . . . . (42) 
n J dq 

by requiring the single-valuedness of the wave function. The integration contour is chosen 
to enclose the two turning points at energy E. Defining the function 
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this quantization condition may be written as 

N SC (E) =n, n = 0,1,2,.... (44) 

The semiclassical energy eigenvalues E sc (n) (n = 0, 1, 2, . . .) are then obtained by inverting 
Eq. Now, the semiclassical expansion of Eq. ([43]) give si 

N SC (E) = \UE) ~\- hJ 2 (E) + 0(h 3 ) (45) 



where 



dq 

Here Jq{E) is the action variable with the property 

dJ (E) 1 



j fjp\ 1 La dS 2k{E,g) 

ME) = ^fd q —-—. (46) 



(47) 



dE u(E) 

Truncation of the series ( |45| ) at order h reproduces the Bohr- Sommerf eld quantization con- 
dition. Note that it is necessary to distinguish the semiclassical quantity N SC (E) from the 
exact N(E) which is a staircase function. For potentials where the series expansion ( f45|) 
terminates at some finite order, N SC (E) smoothly interpolates between the sharp steps of 
N(E). Similarly, the exact quantum mechanical density p(E) has to be distinguished from 
the smooth semiclassical level density 

pUE y.^J^ = ] J--nj i{ E )+ 0^. (48) 

Accordingly, we define a generalized oscillation frequency containing quantum corrections 
by 

Q(E) := j-^ = lu(E) + h 2 J' 2 {E)uj{Ef + 0(h 4 ). (49) 
Another quantity useful below is the energy difference between nearby levels. Defining 
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hAi(E) :=E sc (N sc (E) + l)-E, 
we find with the help of Eq. 



/,,!,( E) = hA { p{E) + '-Af\E) + ^(E) + '^{E) + 0(h 5 ) 



n 2 



I (2) . 



I (3), 



ft 4 



I (4), 



where 



4" 


= Iuj, 








a« 


I 2 

~ 2 


u 2 


■> 






= 6lu 2 4 


/ 3 

+ — LU 
9 


" 21" 


A<?> 


= 121 


2 [uj 3 4]' + 


/4 / r ^ 



Here we have suppressed the explicit energy dependence for clarity. 



(50) 



(51) 



(52) 



B. Semiclassical matrix elements 

As a first step to a calculation of semiclassical matrix elements the normalization factor 
C(E) needs to be determined. Extending previous work on semiclassical matrix elements^, 
one finds to second order in h 



p{E,q) 



l - h 



2 1 dS 2 (E,q) 



u{E) 



+ h 2 C (2 \E) 



(53) 



p(E, q) dq 

where the last equality defines C^ 2 '(E). Below we will need semiclassical expressions for the 
dipole matrix element 

Ql(E) := (E + hAi(E)\q\E). (54) 
It is, however, more convenient to evaluate the momentum matrix element 

P,(£):=(£ + ^(£)N£> (55) 
from which the dipole matrix element can be obtained by means of the relation 
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P l (E) = iMA l (E)Q l (E) 



(56) 



which is a consequence of [H,q] = hp/iM. Combining Eqs. (0), (|55|), fl5"SP, and ( ^T|) one 
finds along these lines 



27riA l (E) 



dq 



' P (E, q) + l 1 ami) _ % .SS 2( E, q) 



2ip(E,q) dq 



dq 



(57) 



x 



exp (i ?) - S (E + hAi(E), q) - h 2 [S 2 (E, q) - S 2 (E + hAi(E), q)}} 



>C(E)p(E, q)C(E + HAi(E))p(E + HAi(E), q) 
A systematic expansion in terms of Ti then leads to 



Qi{E) = Qf\E) + hQ?\E) + \h 2 Qf\E) + h'Qf \E) + 0(h 4 ) (58) 



where the zeroth order term is given by 

Q\ 0) ( E ) = ^hJ d( l ex P HME)t(E, q)} . 



(59) 



Here we have introduced the time 



t(E,q) 



dS (E,q) 
OE 



(60) 



spent by a particle of energy E on its way from the turning point to q. Transforming the 
integration variable and integrating by parts, it is seen that 



Zir Jo 



(61) 



is just the Fourier coefficient of the classical trajectory introduced in Eq. fllTf). The first 
order correction is completely determined by the zeroth order term by means of 

Iuj(E) 



Q { i\E) 



(62) 



and the second order term is given by 
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Q?\E) 



l 2 Uj(E) r 



uQfX{E)- l ^nE)Q?\E) 
Mi 



J' 2 {E) + C^\E)\u{E)Qf\E) 



u(E) 



2n jd q e-^ E ^[2^S 2 {E A ) 
-2u{E)t{E,q)J' 2 {E)- 



2i 



d 



lu(E)p(E,q) dq 
liluj{E)M 2 



S 2 (E,q) 



12p(E,qY 



I 2 ( d 2 



t(E,q)+ uj(E)lu"(E) — 2lo'{E) 2 t(E,q) 



(63) 



In deriving Eq. (|63| ) we have made repeated use of partial integrations. Note that the third 
order coefficient Q\ (E), which is also needed in the sequel, is completely determined by 
Q\ 0) (E) and Qf\E). As a consequence of the symmetry relation 



(n\q\m) = (m\q\n), 



(64) 



or, equivalently, 



Qi(E) = Q-i{E + hAi(E)), 



one finds with the help of Eqs. (|5l|) and (£|5] 



Qi 



(3) 



31 



LO 



Q) 



(2) 



— 31loJ 2 



Q 



(0) 



+ 



LO 



' " (0) 



Qi 



4 



LO 



Q) 



(0) 



(65) 



(66) 



This way we have expressed the quantum corrections to the dipole matrix element up to 
order h 3 by purely classical quantities. 

For the semiclassical expansion of the master equation we need the squared matrix 
element 



B t (E) := (E + hME)\q\Ey 



(67) 



From the above results, its semiclassical expansion is given by 



B l {E) = B?\E) + hB^(E) + '^Bl 2 \E) + ^B^(E) + 0(h^ 



* 2 h " 



where 



(68) 
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(0)1 



B< 2 > 



-u 



2 
I 2 



(0)1 



-UJ 



Q 



(0) 



+ 2Qi 0) Q?\ 



B, 



(3) 



-UJ 



UJ 



+ 3/u; 2 J? 



B, 



(0) 



/ Z 3 



B, 



(0) 



3/ „(2)T 



H 

2 



where we have again suppressed the explicit energy dependence. 



C. Semiclassical expansion of the probability 



(69) 



Let us now assume that there is a continuous function p sc {x) which interpolates between 
the discrete values of p n such that 



Psc(n) = p n , for n integer, 



(70) 



as it is the case for the stationary solution (|30|) of Eq. (^). The probability ([32] ) can now 
be written as 



P(E) = J2 5 ( E ~ E sc (n))p sc (n) = P SC (E) £ 5(N BC (E) - n), 

n=0 n=0 



where we have introduced 



(71) 



P SC {E) = p sc (E)p sc {N sc {E)), 



(72) 



and have made use of Eq. fl45|) and the fact that 



6{E - E sc (n)) = Psc (E)6(N sc (E) - n), 



(73) 



which is a version of the well-known relation 



6(x) = \f'(x )\5(f(x)) 



(74) 



holding for functions f(x) with a single zero Xq. Now, the semiclassical expansion of P(E) 
can be obtained in the following way. With the help of the identity 
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oo oo 

a 



S(x - n) = 8(x) + ^ 8(x - n) (75) 



n=0 n=0 

we can write formally 

oo 1 _d_ -r 

£*(*-") = . + (76) 

n=0 Sinh (3^) 2 

Here the rhs is merely an abbreviation for the series 

oo 9 _ r>2n i 7 

= E ^4r^ 2n = 1 - ^ + ^n x4 + ■ ■ ■ ( 7? ) 

smh(x) ^ (2nj! 6 360 

where the B n are Bernoulli numbers. Now, in view of 

d dE d 



dN, r dN, r dE 



we have 



i a 



-4§r = f 5(N SC (E) - n) = . , 2 ;- ( ?^ 9 . 9(iV sc (E) + i) 



(78) 



2 ^ e(iV sc (£) + -)■ (79) 



" sinh(M-)^) ^ ^ ' 2 
From Eq. fl4"5| ) we find 

9(iV sc (£) + i) = 6(E) - n 2 w(0) J 2 (0)«J(S) + 0(ft 3 ), (80) 

where we have used the fact that Jo(0) = 0. This combines with Eqs. (|79|) and (^8[) to yield 
the semiclassical expansion 



P( E ) n 



2 



Q{E) - —u(E)u(0)5'(E) - h 2 u(0)J 2 (0)5(E) + 0{h A ). (81) 



Psc(E) v ' 24 
Hence, from Eq. ([H]), the probability reads in the semiclassical limit 



h 2 



Q(E) - —u(E)u(0)6'(E) - h 2 u(0)J 2 (0)8(E) + 0(h 4 



(82) 



P(E) = P SC (E) 
Now, as a consequence of Eq. (|30|), we have 

^s ( c q) (^) = ^Psc(^). (83) 
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The equilibrium distribution function then takes the form 



-PE 



(e(E) 



-ti 2 



+ o(h 4 



hZ \u(E) 

Thus, to second order in h equilibrium expectation values (|33| ) are given by 



dEe 



ti' 



-PE 



h 2 J' 2 (E) 



f(E) 



+ ^(0) [/'(0) - 0/(0)] - ftV a (0)/(0) + 0(h*) 



Using this for f(E) = 1, the partition function Z is found to read 

1 



hz = hz r ] — h 2 



/ dEJ 2 (E)e^ E + -MO) + J a (0) 



where 



/Wo 

is the classical partition function. Of course, truncation of the expansion in terms 
the range of temperatures to k B T S> hou as is readily seen from Eqs. (|8^)-d86|). 



D. Semiclassical expansion of the transition rate 



To expand the transition rate in the semiclassical limit we define 



W t (E) := W Nsc{E)+l 



N ac (E) 



which reads in view of Eqs. (^) and 



Wi{E) = —^Bi(E)D{hAi(E)). 
h 



The semiclassical expansion is given by 



wm = ^<\E) + \wf\E) + W?\E) + hW?\E) 



Introducing the cosine moment of the damping kernel 
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7c0) = Rej(iz), 



(91) 



the coefficients of the expansion (|9~0"D take the form 



wl 



(0) 



W, 



(i) 



(2) 



(3) 



2M 




M 

12/3 
+ 



^ + 



B 



(0) 



7c +i/ic'5; o) 7 : 



v 



12B\ 2) ) 7 , 



Z 2 {-6/W + 2/3V}B, (0) -6/3Z 

4cj V + 4^') 2 - 6/3c;V] + 2A\1\uj 2 J' 2 ) B\ Q) + 6\l\ 



|3 u;V 



B 



(0) 



Yc (92) 



+3lW(u') 2 Biy: 



IMu 
24/F 



{l 2 [-4/3 (a/) 2 - 4/W' + 4/3 2 cjo/] - 24/3^} 5< 



?(o) 



+ [l 2 \-2uu" - 2(cj') 2 - 6/W + 2/3V] + 24o;J^} 



B 



(o) 



-6/W 

|3 



(o) 



- 2/ V 



0(0)1'" 



\2f3B\ 2) + 12 



5 



(2) 



7c 



+ 



+ 



2wV" + 2{uo'f + 8cW" - 10/3W) 2 - 4/3wV + 2/3 2 cu 2 ^ 



72o;a;'J2 + 24u 2 J^ - 2A(3u 2 J 2 



.2 tI 



B 



(0) 



+ {|/| 3 [4cu(c/) 2 + 4tA/' - 6/3u;V] + 24|/|cu 2 ^} [b^]' + 12|Z|aABf 2) W c 



+ NZ 4 [Auj(uj'Y + 4wVu/' - 3/3: 
+|Z|V(u/) 3 if>7 c // 



+ 2AlWu'J' 2 } b\ q) + 3/VV) 2 [Sj 



?(0) 



where ^(lu(E)) was abbreviated by 7W and the energy dependence was again suppressed. 



E. Semiclassical expansion of the master equation 



We are now prepared to address the expansion of the master equation (|37|), which may 
be rewritten as 



p (E) = EE i S ( E ~ E n- hME n )) - 5{E - E n )]Wi(E n )p n . (93) 

n=0 l=—n 



Using the relations 
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[5(E-E n -hA l (E n ))-5(E-E n )]f(E n ) = f:(J E ] j [ 5(E - E n )f(E) (94) 



and 

oo oo oo ( co / — 1 



n=0 I=-n Z=l kn=0 n=0 



E E /(«. 0=EiE [/(n, + /(n, -0] - E /(". -0 f > (95) 

n=0 - 

we find from Eq. 

CO / f) \ ( f)\k CO 

W = E^j ^T,{iMm k w l (E) + [A^E)} k w^(E)}p(E) + p B (E) (96) 



where 



P S (£):=-E(^j ^E[^(^)]^(^)E^-^K (97) 
is a boundary term which vanishes since 

W_i(S n )ocB_j(£7 n ) = (n-Z|g|n) = for n < Z - 1. (98) 

Despite the fact that the expansion fl9Tj|) of starts with a term of order h~ 2 , the 

series expansion in h of the rhs of Eq. (|96| ) does not contain terms of order TT l since they 
vanish due to the symmetry of ^(lu^E)) and Q < f > {E) in Z. In leading order we then have 

P(E) = C P(E) (99) 

where 

with A' ' (22) defined in Eq. (|T9|). Hence, we have recovered the Zwanzig energy diffusion 
equation (|18|). 

The next order correction of order h vanishes by symmetry. Thus, the leading quantum 
corrections are of order h 2 . Inserting the expansions provided in the previous subsections, 
one obtains a large number of terms that, after some manipulations, may be combined to 
give 



P(E) = (C + h 2 C 2 )P(E) 



(101) 



where 



2 (3dE 



A^(E) + ^E) + ^ X (E)e^ E ^ E 



d pE u}(E) 
dE 2tt 



+C Q J' 2 {E)u{E). 



(102) 



We also have introduced the functions 



X 



1=1 
OO 



A< 2 > = 2nMu {l 2 [2uJ$ 0) + B\ 2) ] 1c + 2l 3 u 2 J' 2 Bf' ^} 



and 



ttMuj 



1=1 



(103) 
(104) 



1=1 



,3V 



' J 1 Q/6, ,2/, ,/\2 p(0) // 



(105) 



Note that the function ( |104| ) combines with Eq. ( |T9|) to the generalized energy relaxation 
coefficient 



A(E) =AixMn(E)Y J l 2 Bi(E) lc (m(E)) = A^ (E) + h 2 A^ (E) + 0(h 3 ). (106) 



IV. GENERALIZED ENERGY DIFFUSION EQUATION 

The results of the previous section can be summarized as follows. The generalized energy 
diffusion equation containing the leading quantum corrections can be written in the form 

1 d 



P ( E ) = -^{e~ l3E 
K J (3dE\ 



A(E) + h <&(E) 



dE' 







dE" \dE" 2tt J 



P{E) = CP{E) 



(107) 



Here Q(E) is a generalized oscillation frequency containing quantum corrections as defined 
in Eq. fl49|) . A(E) is given in Eq. (|106|) and contains the obvious quantum corrections to 
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the classical energy relaxation coefficient. The functions <f>(E) and x{E) were introduced in 
Eqs. ( |105| ) and ( |103| ), respectively. Apart from the classical quantities 00(E) , Q\ \E) 7 and 
'j(z) that are ingredients of the classical Zwanzig equation, one merely needs to calculate 
Jz\E) given in Eq. ( (4"6| ) and the second order quantum correction Q\ (E) of the dipole matrix 
element given in Eq. ( |B3| ) in order to determine the generalized Fokker-Planck operator in 
Eq. fllOTD explicitly. Since one can show that the equilibrium solution fl84|) is a stationary 
solution of Eq. (|107|) , the only non-trivial quantum correction that survives in equilibrium 
is seen to be the quantity J2(0). 

In the case of an Ohmic bath, i.e. j c (u) = 7, the coefficients become particularly simple. 
Given the functions u>, J2, and B^ 2 \ all quantities entering Eq. (|107|) can be expressed in 
terms of three sums 



A (i) = 4ttM 7 Y, FBI 



2 D W 



Then 



00 



4 H (0) 



4ttM 7 ^ Z 4 ^ 



= u + h 2 oj 2 J2i 



0,2) 



(10£ 



12 y 
lu r 



(u; 2 )> - 2uj 2 <P" 
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A = wA^ + ft 2 



(109) 



1 



1 2 



A. Harmonic oscillator 

The simplest model for a molecular system, which already exhibits many of the charac- 
teristic features of VER, is the harmonic oscillator potential: 

V(q) = \Mujlq 2 . (110) 



For the potential (|110|) one obtains 
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lu(E) = UJq, 
ME) = 0, 



111) 



Qtl(E) 



E 



and the squared matrix elements read 



B ±1 (E) = -L^(E±h^). 
±u ; 2Muol V 2 J 



(112) 



Special to the case of the harmonic oscillator is the fact that Bi(E) vanishes for I ^ ±1. 
The transition rates are 



W ±1 (E) = ±^n(±hu ) (e± 



(113) 



huj y \ 2 

where n(E) was introduced in Eq. (|28|). The Pauli master equation reads in terms of the 
probability P(E) 



P{E) 



7c 



n(huJo) (t 



i\ e + ^ 



-n( 



-hw ) (e^ ^ - l) [ E 



P{E). 



;ii4) 



For the energy relaxation coefficient (|T^) and the coefficients ( |1 03| ) - ( [3~05| ) one finds 



X{E) 



JcE 
j c uj E 
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A i2) (E) = 0, 
$(£) = 0, 



(115) 



so that the generalized energy diffusion equation ( |107|) takes the form 



12 



[(5 2 P'{E) + 2(3P"{E) + P^(E)]E + /3P'(£) + P"(E) 



This equation can be solved in terms of the Fourier transform 



(116) 
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P(u) 



dEe iuE P{E). 



117) 



Despite the fact that P(E) = for E < 0, the integral runs over all energies to collect 
possible contributions of delta functions at E = 0. The Fourier transform of Eq. ( |116| ) takes 
the form 



P(u) 



^wjp(w) + i(/3 - iu)P'(u) 
+—r7r-(l3 - iu) 



12 



(f3 - 2iu)P(u) + u(f3 - iu)P'(u) 



118) 



This first order partial differential equation can be solved by the method of characteristics. 
The general solution is given by 

(3 ~ ( -if3e-^ l u \ 



P(u) 



f3-i(l- e->*)w Po 1/3 - i(l - e"r- t )u / 



+ih 



,(5ujI{1 - e-^)u((3 - iu) \ -i(3 2 ((3 - iu)e-^ l u ~, ( -i^e'^u 



12[p - i(l - e-^)uf \ (3 - i(l - e-^)u \(3 - i(l - e~^ t )u / 
1 r ts o . „,„ + s _ „oi ^ ( —if3e~ lct u \\ 



(119) 



/3-i(l -e-> 4 ), 



[-(1 - e~ 7c > 2 - i(3(3 + e' lct )u + 2/? 2 J P 
where Pq(u) is the initial distribution. For t — > oo this solution tends to 



/3 /5 2 



(120) 



/3 - 24 V (3-iu) 

where we have used Pq{0) = 1 which is just the normalization condition. The inverse Fourier 
transform of this stationary probability gives 



P(E) = e{E)(3 1 + 



tfu 2 f3 2 \ _ 0E % 2 ujI(5 
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24 



((35(E) + 5'(E)) 



121) 



This result coincides with Eq. (|84] ) if one inserts the quantum mechanical partition function 
Z = l/[2sinh(ftu;o/3/2)] and expands up to second order in h. 

As a specific example consider a system with given initial energy Eq, i.e. 



P (E) = 8(E - E ). 



Then moments are readily calculated from Eq. (|119|) using 



(122) 
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u 



P( u ) = l + iU (E(t))-^(E 2 (t)} + ... 



(123) 



This way the expectation value is found to read 



(B(()> = ^ + (i + ?^)( 1 _ 



and the second moment obeys 

1 



(E 2 (t)) 



2 



(f3 2 E 2 - 4(3E + 2)e~ 27c ' + 4((3E - l)e-^ 1 + 2 



2, ,2 



+ 



12 



-4/3£ + 7)e~ 27c ' + 4(/?£ - 2)e~ 7ci + 1 



(124) 



(125) 



From Eq. (|124 ) one finds that for the case of linear dissipation considered here the energy 
relaxation time 



Ti 



1 

7c 



(126) 



is not modified by quantum effects in accordance with the findings by Bader and BerneB. 



B. Morse oscillator 



The fact that the relaxation time is unaffected by quantum effects is of course a special 
feature of the simple harmonic potential (|110| ). As an example of a more realistic model 
for VER we consider a damped Morse oscillator. This system plays an important role in 
the modeling of, e.g., diatomic molecular systems. The Hamiltonian of the Morse oscillator 
reads 



p2 

H =2M + D ° 



I _ e -n(r-r ) 



(127) 



where D is the dissociation energy, r the nuclear distance at equilibrium, and k 1 the range 
of the molecular potential. Here, the P = Mdr/dr is the classical momentum. In terms of 



the scaled variables q = n(r — ro), p = P/ \ / 2DoM, t = tk^2Dq/ M the Hamiltonian takes 
the form 
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where 



H(p,q) = ^- + V(q) 



V{q) = \{l-e-if. 



(128) 



(129) 



Note that now h is dimensionless and stands for hn/^2D M. The classical trajectory for a 
bound particle with energy E < 1/2 is found by solving 



2E - 1 = q 2 + e 



-2q 



2e 



(130) 



with the result 



q(E,t) =ln 



1 + \p2E~ cos (ty/l-2E) 



;i3i) 



1 - 2£ 

Here we have assumed that the particle is at a turning point at t — 0. From this result it is 
seen that the energy dependent frequency is given by 



10 



(E) = y/l - 2E. 



(132) 



To find the Fourier coefficients we write the trajectory in the form 



q(E,t) =ln 



1 + VI - 2E 
2(1 - 2E) 



+ ln[l - z 2 {E)e iu){E)t ] + ln[l - z 2 (E)e 



(133) 



where 



zi /2 (E) 



1 ± VI - 2E 
V2E ' 



(134) 



An expansion of the logarithms gives for I ^ the Fourier coefficients 



Qi{E) 



[ME)} 
\l\ 



\i\ 



.l)V\ 



l-u(E) 
l + u(E) 



i I'l/ 2 



(135) 



Further, the action variable Jq(E) is obtained by integration of Eq. 

rE dE' 



as 



Jo(E) 



1 - Vl - 2E, 



Jo u(E') 

and the second eikonal S 2 (E, q) can also be integrated exactly yielding 



(136) 
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S 2 (E,q) 



-QE + 8E 2 + 6e~ 2q E + 1 - 3e" 9 + 3e~ 2<1 - e~ 3q 



Therefore, according to Eq. 



48E(2E - 1 + 2e~i - e" 2 ?) 3 / 2 

we have J 2 (E) = 0. Similarly, one finds 



C {2 \E) = 0. 



(137) 



(138) 



The contour integral in Eq. (|63|) may be evaluated conveniently by means of the substitutions 



e icu(E)t 

dt 

p(E,t) 
exp[-q(E,t)} 



z, 



dz 



izou(E) ' 

i{z 2 -l){z 2 -l)z 2 
{zz 2 -l){z 2 - z){z 2 2 + l)' 

(4 - 1) 2 ^ 



(139) 



One then finds 

Qf\E) = Qf\E)\l\ 



(zz 2 -l){z 2 - z)(z 2 + l)' 
3lu 2 (E)\1\ 3 - 8lu 3 (E)1 2 - \2\l\E 2 + 2uj\E) 



96E 2 u 2 (E) 



which gives for the squared matrix element 



Bf\E) = { QT\ 2 { E)^m ' + mi^(E) - 2E 2 



(140) 



;i4i) 



12\l\u 2 (E)E 2 

Restricting ourselves to the Ohmic case 7 c (u;) = 7, one only needs to calculate the functions 
( |108|) which read 

A (0) (£) = 2tt 7 [(1 -2E)- 1 ' 2 - l] , 

A (2) (£) = 2n 1 (E + 2)(1 - 2E)- 3/2 , (142) 
(f)(E) = 2^E(l - 2E)~ 3/2 . 

(143) 



Using Eq. ( |109| ) we then find 



X (E) = 2tt 7 ^. 



$(E) = -2tt7 
A {2 \E) = 2tt 7 



E + 2 



4(1 -2E) 2 ' 
E + 2 



(144) 



4(l-2£) s 



25 



Thus, the generalized energy diffusion equation takes the form 



HE) 



h 2 d 



■PE 



_d_ 

dE 



e pE u(E)P(E). 



(145) 



12 BE dE 

We recall that all quantities are dimensionless, in particular the dimensionless energy E is 
measured in units of 2D Q and (3 = 2D Q /kBT. 

Let us now investigate how an initial probability 



P (E) = 5(E - E 



(146) 



evolves for short times. Using 



d 



5 </<s<.)» 



t=0 



dECP (E)f(E) 



(147) 



one finds for the initial change of the energy expectation value and the second moment 




Of 



(E(t)) 



d_ 

Of 



(E 2 (t)) 



t=o 



t=o 



-7 
-27 



a; Jo - (^±^ + lRV(2W - 3)' 



(148) 



It is useful to discuss this result in terms of the relaxation rate R(E ) defined by 

= —R(E a )E . (149) 



t=0 



From Eq. (|148|) one sees that quantum corrections lower the relaxation rate since the mean 
thermal energy is raised by zero point motion. Usually this effect is accounted for by the 
definition of the relaxation time T\ according to 



d 



(E(t)) 



1 



l(E(t)) - (E) 



cqj 



(150) 



Contrary to the harmonic case ujJq ^ E, and the energy relaxation is non-exponential for 
high energies where anharmonicities are important. The usual concept of a relaxation time 
is thus no longer strictly valid. Eq. ( |150| ) may, however, serve in connection with Eq. (|148|) 
to define an initial relaxation rate 1/T±. For the classical relaxation time one then finds for 
(3 3> 1 (that is the potential is almost harmonic for thermal energies) 
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from where it is readily seen that the classical relaxation time increases as the initial energy 
increases, consistent with the fact that the period increases with energy. The semiclassical 
corrections to the relaxation time are then found by using Eqs. fl85| ) and The leading 
order correction to Eq. ( |151| ) is found to read 

Tf - 1 - (I " I) (152) 



Tf ~ ' 12 ' \E J„ 

from where the classical and harmonic limits are readily recovered by setting h = and 
Jo(E Q ) = E , respectively. With the rates defined in this way, we now find an acceleration 
of the relaxation process by quantum fluctuations. 

To investigate the quality of the semiclassical energy diffusion equation in some detail, 
we compare with results of the full master equation. Since for the Morse potential exact 
dipole matrix elements are available^, the initial energy relaxation rate may be calculated 
using 

d °° 

— {E{t))=J2P™(t)W mn {E n -E m ) (153) 

°^ n,m=0 



which is a consequence of the master equation fl22). To compare this with the approximation 



( |148|) , one has to consider an initial distribution 

Pn{0) = <W (154) 

For Ohmic damping we then have 



where E = E. 



noi 



= -R(n )E = 2 7c ^ Eq) _ (155) 

t=0 l=~n c 



ME) = Ilu(E) - l \ (156) 



and B[{E) is given by0 
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M£) _Ar (I , 1 , <AE)\ V (l 4- I - 

Figure [I| compares the relaxation rates R(no) for the parameters h = 0.0943, (3 = 20 (ten 
bound states) and h = 0.0198, (3 = 100 (50 bound states). Since N sc ((E) eq ) = 0.27 and 0.17, 
respectively, in equilibrium only the ground state is populated substantially. Hence, R(n) is 
negative for n = only. The figure shows that for higher states the classical rates agree well 
with the quantum results. For low-lying states, however, while the classical approximation 
fails the semiclassical rates are still in good agreement with the master equation prediction. 

V. CONCLUSIONS 

We have given a systematic semiclassical theory of vibrational energy relaxation based on 
the quantum master equation. Specifically, we have studied a model leading to a generalized 
Langevin equation in the classical limit, which is the basis of much of the earlier work. An 
expansion in powers of H (more specifically, in powers of hido/hsT, where Tlujq/2 is the 
vibrational zero point energy) was shown to yield the Zwanzig energy diffusion equation and 
quantum corrections thereof. For a harmonic potential, the energy relaxation time is not 
affected by quantum fluctuations while it can be reduced for nonlinear potentials. For the 
Morse potential the improved energy diffusion equation was shown to incorporate the main 
quantum effects. 
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FIG. 1. Comparison of energy relaxation rates R(n) for various initial states n based on classi- 
cal, semiclassical and quantum theories for the Morse potential and two sets of parameters. R{n) 
is given in units of j c . Ground state rates are negative since in equilibrium also higher vibrational 
levels are populated. 
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